import numpy as np
from matplotlib import pyplot as plt

r1 = 1.496e8
r2 = 2.279e8
R1 = 6371
M = 1.989e30
M1 = 5.965e24
T2 = 686.98
G = 6.67e-20
a_min = (r1 + r2) / 2
a_max = 2.773e9

a = np.linspace(a_min, a_max, num=100)

p1 = np.sqrt(a) / (np.sqrt(G * M))
p2 = np.pi * a / 2
p3 = a * np.arcsin((a - r2) / (a - r1))
p4 = np.sqrt(2 * a * r2 - 2 * a * r1 + r1 ** 2 - r2 ** 2)

t = p1 * (p2 - p3 - p4) / (60 * 60 * 24)

inside_arctan = np.sqrt((2 * a - r1 - r2) * (r2 - r1) * (2 * a * r1 - r1 ** 2)) / (
    a * r2 - 2 * a * r1 + r1 ** 2
)

act = np.array([])
for item in inside_arctan:
    if item < 0:
        act = np.append(act, np.arctan(item) / np.pi * 180 + 180)
    else:
        act = np.append(act, np.arctan(item) / np.pi * 180)

theta1 = act + t / T2 * 360
theta2 = -act + (2 * np.pi * np.sqrt(a ** 3 / (G * M)) / (60 * 60 * 24) - t) / T2 * 360

plt.subplot(1, 3, 1)
plt.title(r"$t(a)$")
plt.xlabel(r"$a(km)$")
plt.ylabel(r"$t(day)$")
plt.plot(a, t)
plt.plot(a[-1], t[-1], "o", color="black")
plt.grid()
plt.text(2e9, 75, r"$(2.773×10^9, 72.03)$")

plt.subplot(1, 3, 2)
plt.title(r"$θ_1(a)$")
plt.xlabel(r"$a(km)$")
plt.ylabel(r"$θ_1(°)$")
plt.plot(a, theta1)
plt.grid()

plt.subplot(1, 3, 3)
plt.title(r"$θ_2(a)$")
plt.xlabel(r"$a(km)$")
plt.ylabel(r"$θ_2(°)$")
plt.plot(a[:7], theta2[:7])
plt.plot(a[0], theta2[0], "o", color="black")
plt.plot(3.36e8, 500, "o", color="black")
plt.grid()
plt.vlines(a[0], 100, theta2[0], color="black", linestyles="--")
plt.text(a[0], 125, r"$(1.89×10^8, 135.66)$")
plt.vlines(3.36e8, 100, 495.66, color="black", linestyles="--")
plt.text(2.7e8, 505, r"$(3.36×10^8, 495.66)$")

plt.show()
